using System;
using dnAnalytics.Nli;

/// <summary>
/// Matrix vector multiplication. y = alpha*a*x + beta * y, where a is a matrix, and x and y are vectors.
/// </summary>
internal class Program {
	[STAThread]
	private static void Main() {
		//use the MKL BLAS provider
		IBlas blas = Blas.Provider(BlasProvider.Mkl);

		//row ordering
		double[] a = new double[] {1, 1, 1, 2, 2, 2, 3, 3, 3};
		double[] x = new double[] {1, 1, 1};
		double[] y = new double[] {0, 0, 0};
		double alpha = 1;
		double beta = 0;
		int lda = 3;

		blas.Dgemv(BlasOrderType.Row, BlasTransType.NoTrans, 3, 3, alpha, a, lda, x, 1, beta, y, 1);
		foreach(double value in y) {
			Console.WriteLine(value);
		}
		Console.WriteLine();
		
		//this is equivlent to call above, since we are setting the ordering to column and then
		//transposing the matrix.
		blas.Dgemv(BlasOrderType.Column, BlasTransType.Trans, 3, 3, alpha, a, lda, x, 1, beta, y, 1 );
		foreach(double value in y) {
			Console.WriteLine(value);
		}
		Console.WriteLine();

		//column ordering
		a = new double[] {1, 2, 3, 1, 2, 3, 1, 2, 3};
		
		//this is equivlent to call the calls above
		blas.Dgemv(BlasOrderType.Column, BlasTransType.NoTrans, 3, 3, alpha, a, lda, x, 1, beta, y, 1 );
		foreach(double value in y) {
			Console.WriteLine(value);
		}
		Console.WriteLine();
		
		//lets add a row and columns of zeros to a
		a = new double[] {1, 2, 3, 0, 1, 2, 3, 0, 1, 2, 3, 0, 0, 0, 0, 0};
		
		//change lda to 4, since we added a new row
		lda = 4;

		//we still get the same result, since we are telling BLAS to only use the top 3x3 part of a
		blas.Dgemv(BlasOrderType.Column, BlasTransType.NoTrans, 3, 3, alpha, a, 4, x, 1, beta, y, 1 );
		foreach(double value in y) {
			Console.WriteLine(value);
		}
		Console.WriteLine();

		
		/* 
		Output:
		3
		6
		9

		3
		6
		9

		3
		6
		9

		3
		6
		9
		*/
	}
}